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ABSTRACT 

LS I +61 303 has been recently detected as a periodic 7-ray source by the 
Major Atmospheric Imaging Cerenkov (MAGIC) telescope. A distinctive orbital 
correlation of the 7-ray emission was found. This work shows that the range 
of uncertainties yet at hand in the orbital elements of the binary system LS I 
+61 303 as well as in the possible assumptions on the stellar wind of the optical 
companion play a non-negligible role in the computation of opacities to high 
energy processes leading to 7-ray predictions. The geometry influence on the 
propagation and escape of 7-ray photons is explored. With this study at hand, 
we analyse the results of a pulsar wind zone model for the production of 7-rays 
and compare it with recent MAGIC observations. 

Subject headings: X-ray binaries (individual LS I +61 303), 7-rays: observations, 
7-rays: theory 

1. Introduction 

LS I +61 303 is currently one of the most studied high energy 7-ray sources (Albert et 
al. 2006, 2008a,b). Discovered to shine at TeV energies by MAGIC, it shares with LS 5039 
(Aharonian et al. 2005) the quality of being the only two known mildly variable (Torres et al. 
2001) 7-ray binaries that are spatially coincident with sources above 100 MeV listed in the 
Third Energetic Gamma- Ray Experiment (EGRET) catalogue (Hartman et al. 1999). Both 
of these sources show low X-ray emission and variability, and no signs of emission lines or 
disk accretion (see the recent works by Sidoli et al. 2006, Chernyakova et al. 2006, Paredes 
et al. 2007). 
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Extended, apparently precessing, radio emitting structures at angular extensions of 
0.01-0.05 arcsec have been reported by Massi et al. (2001, 2004); but this discovery was 
not confirmed by recent observations (Dhawan et. al. 2006, Albert et al. 2008a). In fact, 
Dhawan et al. (2006) presented observations from a July 2006 VLBI campaign in which rapid 
changes are seen in the orientation of what seems to be a cometary tail at periastron. This 
tail is consistent with it being the result of a pulsar wind. No large features or high- velocity 
flows were noted in any of the observing days, which implies at least its non-permanent 
nature. The changes within 3 hours were found to be insignificant, so the velocity can not 
be much over 0.05c. Albert et al. (2008a) confirmed these finding observing in a different 
period about a year later, and showed that the morphology of the radio emission of the 
system is maintained. This emphasized the possibility that a pulsar is the compact object 
companion, since the changing morphology of the radio emission along the orbit would require 
a highly unstable jet, which details are not expected to be reproduced orbit after orbit. It 
is also to be noted that all other Be binaries known to date have neutron-star companions 
(Negueruela 2004). The reason for this may lie in source evolution if the rotation of Be 
stars is achieved during a period of Roche-lobe overflow mass transfer from its initially more 
massive companion. If the companion loses most of its mass in the process, becoming a He 
star with mass of only a few M Q , when exploding as a supernova it can only leave a neutron 
star behind (see, e.g., Gies 2000, Tauris & van den Heuvel 2006 and references therein). 

There has been a recent burst of activity trying to understand and model the high energy 
multi-messenger emission from LS I +61 303, and especially, but not only, assuming that it 
is a pulsar binary (see, e.g., the latest works by Bosch-Ramon et al. 2006, Bednarek 2006, 
Gupta and Bottcher 2006, Neronov and Chernyakova 2007, Dubus et al. 2006a,b, Romero 
et al. 2007, Torres & Halzen 2007, and Zdziarski et al. 2008). The assumptions about 
the orbital elements and stellar wind properties in these works differ and the impact of their 
assumptions have generally not been explored. Zdziarski et al. (2008) consider dumpiness of 
the polar outflow, as a way to explain the X-ray variability found on timescales much shorter 
than the orbital period. Their analysis of the system also finds that the presence of a young 
pulsar is compatible with all observational constraints. Free-free absorption suppresses most 
of the radio emission within the orbit, including the pulsed signal of the rotating neutron 
star. In here, we do not consider dumpiness of the polar wind but rather the impact (onto 
the computation of opacities to high energy processes leading to 7-rays) of even more basic 
assumptions, like the uncertainties in the orbital elements of LS I +61 303 and the main 
average parameters of the stellar wind of the optical companion. With this study at hand, 
we analyse the results of a pulsar wind zone model (PWZ) for the production of 7-rays and 
compare it with the most recent MAGIC observations. 
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Fig. 1. — System geometry for two different orbital solutions of LS I +61 303, based on 
Grundstrom et al. 2007 (left) and Casares et al. 2005 (right). The phases for Inferior 
conjunction (INFC), Superior conjunction (SUPC), periastron, and apastron are marked 
accordingly to these solutions. The measurements of phases for this system start from the 
radio ephemeris (see. e.g., Gregory 2002), which exhibits periodic radio outbursts, thus, 
phase does not correspond to periastron in any of the orbital solutions. The inclination of 
the system is not shown. Each of the figures are roughly in scale. For the massive star, the 
equatorial wind (disc) is shown assuming a disc radius = 7R S . 



2. System uncertainties and assumptions 



2.1. Orbital elements 



The assumed values for the basic binary parameters of LS I +61 303 needed for our 
study are summarised in Table [TJ They give account of the two recent solutions for the 
orbital elements that are available in the literature (Casares et al. 2005, and Grundstrom et 
al. 2007); which are schematically shown in Fig. [TJ Whereas they are consistent with each 
other concerning some of the parameters or the orbital period; others significantly differ. The 
nominal values for the eccentricity, when comparing Casares et al. (2005) and Grundstrom et 
al. (2007), also span a large range, although both values are barely consistent when looking 
at the given errors: they are quoted as e = 0.55 ± 0.05 for the latter and e = 0.72 ± 0.15 
for the former. The orbital solution by Grundstrom et al. (2007) seems to be in reasonable 
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agreement with the expectations from the radio model of Gregory (2002), which places the 
phase of periastron at <p p =0.33 - 0.40. The longitude of periastron (the angle within the 
plane of the orbit which defines the position of periastron with respect to the observer; e.g., 
for w per = 90°, then the periastron is between the star and observer), significantly influences 
the geometry of the system and results in different orientation of the binary with respect to 
the observer. The two quoted values are u p = 57° ± 9 and u p = 21.0° ± 12.7 (Grundstrom 
et al. 2007 and Casares et al. 2005, respectively). Also the difference for the phase of 
periastron, <p p = 0.301 ± 0.011 and P = 0.23 ± 0.02, plays a role in the interpretation of 
observational data. 

Grundstrom et al. (2007) obtained 100 spectra of LS I +61 303, all of them showing the 
Ha emission line and the He I A6678 feature. To obtain the orbital solution, they assume 
that the radial velocity variations of this composite profile represent the motion of the Be 
star, since the line formation probably occurs very close to the photosphere of the Be star 
itself. Unfortunately, even when the number of observational radial velocity data points 
have been much increased by Grundstrom et al. (2007), the paucity of measurements in the 
orbital phase range 0.4 - 0.5, and the possibility that the Balmer lines could be contaminated 
by the disc of the star, still cast some doubts as to what is the real orbital solution of the 
system. 

In what follows, we will work with both sets of orbital solutions and show that such 
differences, and others given in Table [U produce non-negligible changes in the opacities of 
leptonic processes that lead to high energy emission. Important geometrical magnitudes, 
such as the angle to the observer as measured from the pulsar site as a function of phase 
along the orbiio are notably affected by the choice of orbital solutions, since they depend 
directly on the geometry, as can be seen in Fig. [2j In turn, this and other differences 
unavoidably lead to changes in opacities. 

The orbital period of LS I +61 303 is ~ 26.5 days (Taylor & Gregory, 1982), with 
the best measurement coming from the radio observations (26.4960 + 0.0028) days (Gregory 
2002). The massive star is of B0 Ve type (Hutchings & Crampton 1981, Paredes et al. 1994). 
Based on the spectral type (e.g., Cox 2000), its radius and temperature are R s ~ 10i? Q and 
T s = 2.8 x 10 4 K, as used by Casares et al. (2005). If, as in Grundstrom et al. (2007), we 
rather assume a radius of R s ~ 6.7R Q (e.g. Harmanec 1988), the effective temperature is 
T s = 3.0 x 10 4 K. A different temperature changes the photon field against which we have 
computed inverse Compton processes. The dimensions of the star have also an impact on the 
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presented results for each orbital configuration yielding to the change of relative separation, 
what is shown in Table HJ the massive star is not a point-like source, what we have also taken 
into account. 

The estimation of the inclination of the orbital plane is related to the knowledge on the 
mass of each of the stars. To give an example; from the mass function, f(M) = 0.0107 M & 
(Casares et al. 2005), and taking the possible inclination as i = 30° or 60° with the mass of 
the massive companion in the middle of their preferred mass range M s = 12.5 M Q , one can 
get the mass of the compact object within m ~ 2.5 — 1.5 M & (higher inclination implies a 
lower mass of the compact object, favouring a pulsar scenario). Based on these estimations, 
we get the semimajor axis of the binary as a ~ 0.42 AU (a ~ 6 x 10 12 cm). Casares et 
al. (2005) conclude that the compact object would then be a neutron star for inclinations 
25° < i < 60° and a black hole for i < 25°. Finally, a pulsar scenario for the binary implies 
an assumption on the spin down luminosity of the pulsar (thus, a free parameter of the 
model). Typically, for young pulsars it is found in the range 10 36 - 10 37 erg s -1 . We set this 
parameter to L s d = 5 x 10 36 erg s _1 for the following estimations. 



2.2. Stellar wind 

The stellar wind of the Be companion is assumed to have two components, one related 
with a polar contribution, and the other, with the equatorial disk (e.g. Waters et al. 1988). 
The polar wind is radiatively driven (e.g. Castor & Lamers 1979), for which the velocity law 
is V p {r) = VP° l + (V£ o1 - V po1 ) (1 - R s /rf , and where typically, (5 « 1, V po1 = O.OlVg* and 
the terminal velocity of the wind V^ 1 = 1500 — 2000 km s -1 , with mass loss rates obtained 
from UV resonance lines: M p ~ 10 _8 M Q yr~ 1 (e.g., Snow 1981, Waters et al. 1988). For 
the equatorial wind (e.g., see Waters 1986, Waters et al. 1988, Gregory & Neish 2002), the 
velocity law is V eq (r) = V^ q (r/R s ) m where m = 1.25 (from the density profile assumption 
for LS I +61 303, see Waters et al. 1988) and Vq 9 = 5kms _1 . These authors further assume 
that the terminal velocity of the equatorial wind is ~ few 100 km s _1 , so that <C V^ 1 , 
an assumption mimicked by many other works in the field. The velocity laws of both the 
polar and equatorial outflows are depicted in the left panel of Fig. [31 

The mass loss rate for the equatorial wind in Be stars is generally obtained from IR excess 
and Ha lines. For LS I +61 303 it was quoted as M d = (1 —4) x 10~ 7 M Q yr _1 by Waters et al. 
(1988), although uncertainties in this value exist, since it depends on additional knowledge of 
the disc half opening angle and initial wind velocity. A relationship between the mass fluxes 
of the two wind regions is (e.g., Lamers &Waters 1987) Fa/F p w Mm/ (Muv sin 9), where 
F<i/F p is the ratio of the mass fluxes of the equatorial and polar regions, respectively. Using 



- 6- 



Table 1: LS I +61 303: system parameters and a comparison with LS 5039 



Parameter 




Adopted value 


LS 5039 


Distance to the system D 




2.3 kpc 


2.5 kpc 


Semimajor axis a 




6.3 X 10 12 cm 


0.15 AU ~ 3.5 R s 


Mass of star M s 




12.5 Mq 


23 Mq 


Grundstrom et al. 2007 Casares et al. 2005 


Eccentricity of the orbit e 


0.55 


0.72 


0.35 


Longitude of periastron u p 


57° 


21° 


226° 


Phase of periastron tf> p 


0.301 


0.23 





Periastron separation d p 


~ 6R S 


~ 2.6 Rs 


~ 2.3 Rs 


Apastron separation d a 


~ 21 R a 


~ 15.7 R s 


~ 4.7 Rs 


Radius of star R s 


6.7 R Q 


10.0 R e 


9.3 Rq 


Temperature of star T s 


3.0 x 10 4 K 


2.8 x 10 4 K 


3.9 x 10 4 K 


Mass loss rate of star (polar wind) M p 




10" 8 M yr" 1 


lO- 7 M yr- 1 


Wind termination velocity (polar wind) Vj£° 




2000 krns" 1 


2400 km s" 1 


Wind initial velocity (polar wind) Vq° 1 




20 kms -1 


4kms" 1 


Mass loss rate of star (equatorial wind 1) 




1.3 X 10- 7 M Q yr- 1 




Mass loss rate of star (equatorial wind 2) M<j 




1.3 X lO" 6 M yr" 1 




Wind termination velocity (equatorial wind) V^J 




300 km s" 1 




Wind initial velocity (equatorial wind) Vq 




5 km s — 1 




Rotational velocity of the star V ro t 




400 km s" 1 




The half opening angle of the disc 9^ 




15° 




The radius of the disc 




7 R s 






Fig. 2. — The angle to the observer, measured form the pulsar site, as a function of phase 
along the orbit for two different values of the binary inclination angle i The angle was 
calculated for two different orbital projection, based on Grundstrom et al. 2007 (G - left) 
and Casares et al. 2005 (C - right). The phases for inferior conjunction (INFC), superior 
conjunction (SUPC), periastron, and apastron are marked accordingly to the model. The 
angle is defined to be a b s = for propagation outside with respect to the massive star. 
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Fig. 3. — Left: the velocities of the two components of the massive star wind as a function 
of distance from the star centre. The parameters of the wind are in Table [TJ Right: kinetic 
power of the polar and equatorial massive star winds versus radius (from the centre of 
the massive star). For polar wind the power is calculated for mass loss rate M§ = M p = 
10~ 8 M yr -1 and M 9 = M p = lO _9 M yr~ 1 , for comparison. For the equatorial wind 
the power is calculated for the scaling factors fi = 50, = 500 and M 8 for both cases. 
The shaded areas correspond to the range of the binary separation along the orbit (from 
periastron to apastron) for two models of the binary: based on Grundstrom et al. 2007 
(light blue) and Casares et al. 2005 (light yellow). The distances which correspond to the 
separation at the periastron and apastron for each model are marked with vertical blue lines. 
The pulsar spin down power is also marked. 

F d = M IR /(4nr 2 sin 6) and F p = M uv / (47rr 2 ) , with an opening angle of 15°, typical values 
are = (10 — lO 3 )M p sin0. In what follows, for the models investigated here we use the 
description = fM p sin 6, where / is the specific scaling factor defining the model. The 
difference between the mass loss rates of the polar and disc regions should be significant to 
have a luminosity in X-rays as observed in most of Be stars. But in the case of LS I +61 303 
as there is no accretion signatures, this value is difficult to set. In addition, based on results 
from X-ray and IR observations, it was found that the densities of the Be star disc region 
are higher than in the case of isolated Be stars (e.g., Reig et al. 2000). An accompanying 
neutron star can trim the disc, preventing further growth, and therefore making it denser. 
Estimations of the disc radius from Ha lines have been made by Grundstrom et al. (2007), 
finding it in the range of ~ 5R S . In what concerns the high energy 7-ray emission, the 
value of the disc truncation is not a trivial parameter. The smaller the radius at which the 
equatorial wind is terminated, the more time the pulsar spend outside this internal disc. In 
the following, we assume the disc radius ~ 7R S , to account for the fact that the disc can 
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occasionally grow to larger radius (e.g., see Grundstrom et al. 2007, also Zdziarski et al. 
2008). Still, the disc size is much smaller than typical values found for isolated Be stars. The 
pulsar is then subject to the equatorial wind only during part of the orbital period, and we 
consider it subject to the polar component otherwise. In Table |5] we put the characteristic 
phase and angles (measured from the periastron), where notation "disc IN" corresponds to 
the phase when pulsar enter into the equatorial wind region, and "disc OUT" when it leaves 
this region. Intrinsic to the uncertainties in these parameters is the fact that for such wide 
choice of polar /disc features we may have a variety of possible situations: the polar wind 
mass flux can be stronger than the equatorial one, or viceversa. To complicate things, the 
ratio of mass fluxes may change with separation as is due to different velocity laws in the 
equatorial and polar winds. 

The fact that the pulsar is not in the disk region all the time give rise to mixed wind 
models, where the influence of polar and disk wind regions must change along the orbit. To 
explore this effect, we assume that the difference in the mass loss rates is in the range of 
/ ~ 50 — 500 what gives (together with the assumption for half opening angle 9 = 15°) the 
mass loss rates themselves as Md ~ 13 — 130 x M p . These values are summarised in Table [U 
where they are referred to as model 1 (/ = 50) and 2 (/ = 500). In the right panel of Fig. [3] 
we see how the power of the massive star wind changes with respect to the assumed constant 
pulsar power along the orbit. Depending on /, we get two quite different scenarios as the 
power of the equatorial wind for model 1 (/ = 50) is lower then the pulsar spin down power 
only for large distances from the massive star, while for model 2 (/ = 500) it is larger than 
L s d already at a few stellar radius. We remark that these scenarios are all within current 
uncertainties in the knowledge of the system. 

Table 2: LS I +61 303: characteristic orbital phases in different geometries 



phase 


Grundstrom 2007 


Casares et al. 2005 


angle 


°] phase 


angle 


o 


phase 


SUPC 


213° 


0.035 


249° 




0.161 


periastron 


0° 


0.301 


0° 




0.23 


INFC 


33° 


0.324 


69° 




0.257 


apastron 


180° 


0.801 


180° 




0.73 


disc IN 


311° 


0.265 


239° 




0.141 


disc OUT 


49° 


0.337 


121° 




0.319 
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2.3. Hydrodynamic balance 

Based on the hydrodynamic equilibrium of the flows, the geometry of the termination 
shock is described by the parameter 77 = MiVi/ M V , where MiVi and M V are the loss 
mass rates and velocities of the two winds (Girard and Wilson, 1987). If one of the stars 
is a pulsar of a spin down luminosity L S( i and the power of the massive star is M S V S , the 
parameter rj can be calculated from the formula rj = L s d/(c(M s V s )) (e.g., Ball and Kirk 
2000). This hydrodynamical construct only fixes the shock position, and since we study 
PWZ processes in this paper, the current approximation is sufficient: In this approach, the 
shock morphology is not dealt with and plays no significant role in the result. Note that for 
f] < 1, the star wind dominates over the pulsar's and the termination shock wraps around 
it. Note also that for rj = 1, the shock is at equal distance, d/2, between the stars. As a 
first approximation, the shock will be symmetric with respect to the line joining two stars, 
with a shock front at a distance r s from the pulsar r s = D^frj/ (1 + y^). The surface of the 
shock front can be approximated then by a cone-like structure with opening angle given by 
ip = 2.1 (l — ?7 2 / 5 /4) 77 1 / 3 , where fj = mm(f|,f| _1 ). This last expression was achieved under 
the assumption of non-relativistic winds in the simulations of the termination shock structure 
in Girard and Wilson (1987), albeit it is also in agreement with the relativistic winds case 
(e.g., Eichler and Usov, 1993; Bogovalov et al., 2007). 

7/-values for the assumed magnitudes of the LS I +61 303 system are depicted in Fig. 
IH (top panels). The values calculated in the case of a polar wind only are marked by a 
thin dotted line (mostly covered by thicker model lines, since the pulsar is out of the disc 
most of the orbit) E Assuming the parameters of the equatorial wind only, the //-values are 
marked by thin dashed (green) line for model 1 with / = 50, and thin solid (green) line 
if / = 500, these are also plotted beyond their grange of validity to facilitate interpreting 
how the mixed models 1 and 2 are constructed. With the disc radius = 7R S and the 
half-opening angle $d = 15° the phases when the pulsar passes trough the disc region are 
shown as shaded area. We can see that the shock front changes position between the stars 
with the orbital phase (depending on the separation), and that 77 for different models can be 
larger or smaller than 1. In general (for both binary orbital scenarios), the wind interaction 
in the polar wind regions gives rj > 1 meaning that the shock is closer to the massive star 
than to the pulsar. It can be seen that assuming a given power for the equatorial wind (not 
only the pulsar's) is essential to any model. 

Regarding spatial information, we see that for different parameters of the equatorial 
wind, rj changes such as for model 2, the PWZ is terminated closer to the pulsar, see the 



2 The description that follows corresponds to the on-line version of the figure. 
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Fig. 4. — Top panel: The r\ parameter for two models of the LS I +61 303 binary, based 
on Grundstrom et al. 2007 (G - left) and Casares et al. 2005 (C - right), model 1 (/ = 50, 
dashed line), model 2 (/ = 500, solid line). As a comparison the dependencies for the polar 
and two equatorial winds (see Table [I]) are shown (blue dotted line for polar wind and thin 
green dashed and solid lines for equatorial wind in model 1 and 2, respectively ). See the 
text for a full explanation. Bottom panel: The stand-off distance of the shock (measured 
form the massive star centre), r so d, in units of stellar radius R s , as a function of binary phase 
for two models of the system (G, C). 
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bottom panels in Fig. H] which show the stand-off distance corresponding to those scenarios 
(the stand-off distance is independent on the binary inclination as it is measured along the 
separation axis from the massive star centre). In these bottom panels, we mark results 
corresponding to model 1 (/ = 50) with a dashed line and to model 2 (/ = 500) with a 
solid line. The stand-off distance calculated in the case of polar wind only (isotropic wind) 
is marked by thin dotted line (again mostly covered by model lines, since the pulsar spends 
most of the orbit out of the disc). Assuming only the parameters of the equatorial wind the 
stand-off distance is marked by a thin dashed (green) line for the model with / = 50, and 
with a thin solid (green) line if / = 500. INFC, SUPC, periastron, and apastron phases are 
marked. 

Note that for orbital parameters from Casares et al. 2005 (right panels), the separation 
at periastron is very small, and the values of rj are in fact not a good approximation for 
such geometry (the stand-off distance can be less than the distance from the massive star 
surface). In this case, the shock is assumed at a minimal distance from the massive star, set 
at 1.3R S (see the dips around periastron already in the top- right panel of Fig. 0J. This will 
also impact on the distance to the shock from the pulsar side, see Fig. |5j 

For investigating the high energy photon production already in the PWZ, an important 
parameter is the distance from the pulsar to the front shock in the direction to the observer, 
r s h, (see Sierpowska-Bartosik & Torres 2007, 2008b). This magnitude, differently to the 
stand-off distance, depends on the orbital inclination. In Fig. [5] we show this distance along 
the orbit for both scenarios of the binary geometry, two inclinations (see the differences in 
angle to the observer in Fig. [2]), and the two models considered for the stellar wind (different 
line styles in each plot). For comparison, the lines for different wind parameters (polar and 
equatorial) are plotted all along the orbit. The distance to the shock for models 1 and 2 are 
consistent with the values for isotropic polar wind (thin dotted blue lines in each figure) for 
phases when the pulsar is in its corresponding region of influence. 

It is interesting to compare the phase periods when the PWZ is not terminated in the 
direction to the observer (i.e., an electron propagating in this direction will not find the 
shock) based on different models and geometries. In Grundstrom et al. (2007) geometry 
(G) and within model 1, we have long period of a non-terminated PWZ which changes with 
the system inclination: it is ~ 0.22 — 0.73 for i = 30°, whereas it is ~ 0.24 — 0.56 for 
% = 60°, what comes from the interaction of the pulsar wind with the polar wind only (no 
termination shock in the equatorial wind; the efect is only in the direction to observer). For 
the same scenario but within model 2, the equatorial wind is relatively stronger (rj < 1) 
and the pulsar wind terminates in addition to previous phases in the disc region at short 
periods prior to periastron <fi ~ 0.26 — 0.28 (i = 30°) or <p ~ 0.26 — 0.27 (i = 60°). In the 
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Fig. 5. — The distance (in units of stellar radius R s ) from the pulsar to the termination shock 
(r s h) in the direction to the observer as a function of binary phase, for different Be wind 
models (where r d = 7R S and Od = 15°) and geometry based on Grundstrom et al. (2007, 
top, marked with G) and Casares et al. (2005, bottom, marked with C) and two inclination 
angles: i = 30° (left), i = 60° (right). The termination shock is marked with thick dashed 
line for model 1 and thick solid line for model 2 (Tabled]). r s h was also calculated in the case 
of polar and equatorial winds only to facilitate visualising the components of model 1 and 2. 
Assuming the parameters of the polar and equatorial winds as for isotropic wind the shock 
distance is marked by thin blue line for polar wind, dashed green (marked as "equatorial 1") 
line for 1 model with / = 50, and solid green ("equatorial 2") line for / = 500. The results 
from polar wind parameters are covered most of the phases with the lines for mixed models. 
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case of geometry based on Casares et al. (2006) we have shorter termination epochs. For 
i = 30° in model 1 the wind is terminated only for ~ 0.04 — 0.14 (what corresponds to 
the pulsar wind interaction with the polar wind). In model 2 this period becomes longer, 
<p ~ 0.04 — 0.22, due to relatively stronger wind in eqatorial region. When the inclination 
is set to % — 60° the magnitude of the angle to the observer is larger what results in shorter 
periods for the non-terminated pulsar wind. Then for both models, 1 and 2, we have the 
same phases when the wind is non-terminated, <fi ~ 0.22 — 0.88, with differences in the value 
of r s h. Note that in that last cases the radius r s h increases from few stellar radii to infinity 
just prior to the periastron. 



3. Opacities 

The optical depths to inverse Compton interaction and 77 absorption were calculated 
for these models of the binary and its geometry to investigate high energy photon production. 
We used the full Klein-Nishina cross section (see the Appendix of Sierpowska-Bartosik & 
Torres 2008 for full details). We have computed the opacities for particles propagating to the 
observer from the position of the pulsar (assumed in the model to be the place of injection). 
We studied the influence of the different orbital elements (Grundstrom et al. and Casares 
et al.), the same inclinations of the orbit (60°, in agreement with Hutchings & Crampton 
(1981) who actually favoured higher values of i), and assumptions on the stellar winds from 
model 1, see Fig. [61 The results constitute a set of several plots that give account of the 
changes in opacities produced by uncertainties in our knowledge of fundamental parameters 
of the binary. 

Optical depths above unity are only found for low energies of pairs for inverse Compton, 
whereas for both pair production and inverse Compton processes with higher initial energy, 
interactions are less probable. A high inclination in the orbital element solution provided by 
Casares et al. (2005) provide the highest opacities, as the angle to the observer reaches the 
maximum close to periastron (see SUPC position and also Fig. [2]) where the radiation field 
density is the highest. This dependence of the angle to the observer along the orbit, together 
with the fact that in Casares et al. solution the system is more compact around periastron 
(larger radius of the star, larger eccentricity) than in Grundstrom et al.'s, results in higher 
optical depths for phases when the pulsar crosses the disc region. The termination of the 
pulsar wind do not influence the opacities much, as in the geometry of Grundstrom et al. 
and within model 1, the PSR wind is unterminated also in the disc region (<j) ~ 0.24 — 0.56), 
while in Casares et al.'s solution the wind is terminated up to phases close to periastron 
(termination for <p ~ 0.22 — 0.88). We find a similar behavior of the opacity curves between 
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SUPC and INFC in both geometries, so the most important influences can be ordered as 
follows: geometry (angle to the observer + binary separation) - radiation field - and then 
the influence from the shock termination (i.e., the jump at the disc we see in C comes from 
change of the r s h of the order of a few radii) 

Of course, the optical depths depend on energy so that, generally but not always, the 
top curves (larger opacities) correspond to the lower energy of injected particle (in these 
cases, E = 100 GeV). The optical depths for pair production, in particular, present a peak 
in energy, which location depend on the angle of propagation. In this case, then, we may 
have lower optical depths for the lower energies explored- see the phases around periastron. 
In summary, for the chosen particle energies: the optical depths for E = 100 GeV on ICS 
are larger than those of pair production, for E = 1 TeV they are both at a comparable level, 
and for E — 10 TeV the optical depth for pair production are larger. 

Notice that in a scenario where photons are produced by pairs accelerating along the 
shock, where magnetic field influence trajectories, photons are produced at different 
directions, not only at the specific observer's angle, and from different places of the shock. 
This would require a full 3D treatment of the cascading process on which we will report 
elsewhere. But in order to describe the general dependencies of the emission to the densities 
and anisotropics local to the interaction regions, we present here the opacities in the direction 
to the observer at the specific phase and stand-off distance. 

The optical depths at stand-off distance for model 1 are shown in Fig. [7J A more 
complicated scenario appears in the case of Casares et al. (2005) orbital solution (right 
panel of Fig. [7J, where the changes due to a minimal distance of the shock from the massive 
star for a limited-phase range around periastron, are seen. For this model, the optical depths 
for phases within the equatorial wind are higher as the stand-off distance is closer to the 
massive star. In the case of inclination 60°, as we depict, there is a double-peak structure 
in the opacities, because for these limited phase-range, the directions to the observer (when 
starting from the stand-off distance) are eclipsed by the star surface; and in that cases the 
optical depth is lower (dip in the curve) since it is calculated only to the star surface. This 
double peak-structure is entirely geometrical, and is absent in the results for the opacities 
for inclination 30°, as is also in the case of the Grundstrom et al.'s (2007) orbital solutions. 
In this latter case, the influence of the disc is not significant mainly because the orbit occurs 
with a larger separation of the system. 

Only in case of Casares et al. (2005) parameters and high inclination values, the optical 
depths for both processes change significantly in the disc region with respect to the polar 
wind phases (the smallest and biggest angles of propagation happen for relatively small 
separations close to periastron). 
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Fig. 6. — The optical depths for Inverse Compton scattering (black lines) and 77 absorption 
(grey lines) calculated for particle injected at the binary separation distance and propagating 
up to the termination shock in the direction to the observer for different primary energy. 
The optical depths were calculated for the model 1 with / = 50 (r^ = 7, R s and 6d = 15°) 
and inclination i = 60°. The left panel correspond to the orbital solution by Grundstrom et 
al. (2007); the right one, by Casares et al. (2005). The electron/photon primary energy is 
set to 100 GeV (solid lines), 1 TeV (dashed lines) and 10 TeV (dotted lines). 



It is clear that the uncertainties in the system play a non-negligible role in the opacities, 
this in turn will translate into complexities of lightcurve and spectral evolution along the 
orbit that can be tested with future quality of data. For futher investigation we choose the 
geometry based on Grundstrom et al. (2007) and model 1 for the massive star wind, i.e., 
moderate values for the mass loss rates in the equatorial wind. 

Fig. [S] shows the optical depths for 77 absorption and Inverse Compton scattering 
calculated at all places in the orbital plane of the system, for ITeV particles traveling in 
the direction to the observer. The orbit of the pulsar, the star (with its physical size), and 
the position of the shock are marked therein. Fig. [9] shows something more complex yet: it 
presents a comparison of the for ITeV particles maximal optical depths for Inverse Compton 
scattering calculated at all places in the orbital plane of the system (i.e., particles traveling 
in a direction tangent to the physical size of the star), with the optical depths found at the 
position of the shock and the pulsar orbit when particles are traveling in the direction to the 
observer. We see that the differences in optical depth produced by the direction of movement 
are clearly important; and thus the system geometry will decisively impact the lightcurve. 
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Fig. 7. — The optical depths for Inverse Compton scattering (black lines) and 77 absorption 
(grey lines) calculated for particle injected at the stand-off distance and traveling in the 
direction to the observer for different primary energy. The optical depths were calculated 
for the model 1 with / = 50 (r^ = 7, R s and 64 = 15°) and inclination % = 60°. The left 
panel correspond to using the orbital solution by Grundstrom et al. (2007); the right one, 
by Casares et al. (2005). 
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Fig. 8. — Optical depths for 77 absorption and Inverse Compton scattering calculated at all 
places in the orbital plane of the system, for ITeV particles traveling in the direction to the 
observer. The Grundstrom et al. (2007) orbital solution is used. The orbit of the pulsar, 
the star (with its physical size), and the position of the shock are marked. 



12.00 



3.16 




-10 -8 -6 -4 



Fig. 9. — Comparison of the for ITeV particles maximal optical depths for Inverse Compton 
scattering calculated at all places in the orbital plane of the system (i.e., particles traveling 
in a direction tangent to the physical size of the star), with the optical depths found at the 
position of the shock and the pulsar orbit when particles are traveling in the direction to the 
observer. The Grundstrom et al. (2007) orbital solution is used. The orbit of the pulsar, 
the star (with its physical size), and the position of the shock are marked. 
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4. Application to a PWZ model of 7-production 

We apply these results to investigate the scenario for VHE photon production in the 
PWZ. As we have shown in previous papers (Sierpowska-Bartosik & Torres, 2007,2008b) the 
geometry of the system is very important in this scenario because of, among other things, 
the linearity of particles' propagation (and further on cascading) from the injection place. 
As an example for LS I +61 303, the PWZ model for the system was investigated based on 
the parameters of model 1 for geometry scenarios based on both orbital solutions, Casares et 
al. and Grundstrom et al.; in both cases, with an inclination fixed at % = 60°. Together with 
the velocities of the wind in the two regions as given in Table [TJ this scenario is such that 
the interaction of the pulsar and massive star winds set the shock almost at half-separation 
distance (with r\ close to 1), except for a very limited range of phases around the periastron 
where the shock is closer to the massive star, see Fig. HI and consequently, where most of 
the difference between models are found. 

The pairs are assumed to come from within the PWZ, and we consider those in the 
direction to the observer, see FigfSJ The electrons propagate linearly and interact with the 
soft photons from the massive star (anisotropic radiation) via inverse Compton. Photons 
following the same direction of the initial electron are thus created, which can initiate further 
pair production in the same photon field. These cascades are followed up to the termination 
shock, where the electrons are trapped by the local magnetic field but the photons pass 
by. Originated at the PWZ, photons can be absorbed in the massive star wind region 
(MSWR). All these processes are included in the simulations, see Sierpowska-Bartosik & 
Torres (2007,2008, 2008b) for details. In this model, apart from the orbital solution and 
wind parameters, we need knowledge of the pair population interacting with the stellar 
photon field, as well as the normalisation factor (the amount of pulsar power transferred to 
primary pairs). We assume, based on both, observational results from high energy 7-rays 
(Albert et al. 2006, 2008a,b) and kinetic plasma physics studies in conditions found for the 
PWZ of pulsars (e.g., Jaroschek et al. 2008 and references therein) that this population 
is described by a power law. The fraction of the pulsar spin-down power ending in the 
interacting pairs can then be written as: TL S( i = J N e + e -(E)EdE. Assuming that the 
distance to the source is d — 2.0 kpc, the normalisation factor needed to compute the number 
of electrons travelling towards Earth is A = N e + e - /And 2 . In the case of a power-law in energy, 
N e + e -(E) oc E~° H . The photon spectra are simulated for 20 phases chosen along the orbit of 
the binary (avoiding the phases close to the change between the polar and equatorial wind), 
which allows to compute the high energy lightcurve. 

At TeV energies, the MAGIC observational lightcurve presents a broad maximum cor- 
responding to phases from 0.5 to 0.8. Looking at Casares et al. and Grundstrom et al. 
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geometry (Table [2]) we can see that this corresponds to phases prior to apastron (and it is 
not connected with an INFC-SUPC distinction, as in the case of LS 5039). Interestingly, this 
corresponds to phases where the PWZ in the direction to the observer is unterminated, and 
an increasing binary separation with the pulsar already outside the disc region. For the TeV 
maximum broad phases (between 0.5 and 0.8) we assume an initial spectrum of electrons 
given by a slope a = —2.6 and a normalization T = 0.1. The MAGIC data at other phases 
along the orbit are not particularly constraining, since in reality they represent only 95% 
CL upper limits at high levels of flux (Albert et al. 2008b). To explore how this model can 
fit the MAGIC light curve (particularly if we assume very low values of fluxes at periastron) 
we change the parameters of the initial power law spectrum for electrons along the orbit up 
to a = —2 and T = 0.01; implying, within the framework of this model, orbital variations in 
the relativistic population of particles, something already suggested in the case of LS 5039 
(Sierpowska-Bartosik & Torres 2008, Dubus et al. 2008). Needed variations in this case, are 
however stronger than those of LS 5039. 

Results for the two orbital models studied are qualitatively similar, but differ in the 
details, especially in what concerns to the lightcurve. These differences in the details are 
only coming from that produced by the uncertainty in the orbital solution elements, given 
that we have kept the same parameters in the rest of the variables (inclination, electron 
population, stellar wind, etc.). We show these results in Figure [TO]. In there, we also show 
with thin lines the results we would get for the phases of the central VHE maximum with a 
5% increase in normalization. Again, we emphasize that the results of the model at phases 
where there is no reliable detection are unconstrained, and are chosen to show a very low 
level of fluxes (larger fluxes are easy to achieve possible). Further constraints both by Fermi 
and MAGIC II in these orbital phases are essential to fix the properties of the model. 



5. Concluding remarks 

This work shows that the range of uncertainties yet at hand in the orbital elements of 
the binary system LS I +61 303, as well as in the possible basic assumptions on the stellar 
wind of the optical companion, play a non-negligible role in the computation of opacities 
to high energy processes leading to 7-rays. It is not only the face value of opacities what 
changes, but the geometrical influence on the propagation and escape of 7-ray photons is 
also different depending on the assumed parameters and orbital solutions. Detailed models 
of the binary needs to take these differences into account when making precise predictions 
for 7-ray fluxes. Hopefully, new observational studies will help to clarify at least some of 
the variables in the orbital solution. As an example, we have shown that a conceptually 
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Fig. 10. — The gamma ray spectrum (where the TeV maximum occur, cf) ~ 0.65) and 
lightcurve from interacting electrons in the PWZ, for system parameters based on the orbital 
solutions by Casares et al. (2005) - C and Grundstrom et al. (2007) - G. See the text for 
details. MAGIC data points are separated between those having less than 2a confidence (for 
which upper limits are also given) and higher confidence ones (shown as filled circles). Data 
is taken from from Albert et al. (2008b), error bars are statistical only. 

simple PWZ model for the production of 7-ray emission (which details are discussed more at 
length in previous works) is qualitatively close to the data. If admitting orbital variations of 
the interacting electron populations along the orbit this model can also predict low level of 
fluxes for phases where MAGIC observations have up to now provided (yet non-constraining 
upper limits). Although this model is qualitatively in agreement with the MAGIC results, 
it is, however, not complete. Particularly in the case of LS I +61 303, we need to consider 
the acceleration of electrons at the shock front itself: particularly for phases when opacities 
in the PWZ are not so high and many electrons will arrive there even if they are accelerated 
within the PWZ. Fig. [TT] presents a comparison of the corresponding optical depths in the 
LS I +61 303 system leading to the spectrum and lightcurve shown in Fig. [10] with the 
opacities found for the LS 5039 model (Sierpowska-Bartosik & Torres 2008b). To have a 
more direct comparison of the values we show the LS I +61 303 opacities after shifting 
the phase of periastron to <p = 0. With no subsequent shock contribution, the fraction of 
spin-down power ending in relativistic electrons needed to fit the MAGIC TeV maximum 
(10% L sd ) is larger than that found for LS 5039, albeit still possible. Subsequent cascading 
in the MSWZ, if many electrons are reprocessed at the shock, can not be neglected. We 
expect to report on a more detailed model of gamma-production taking into account these 
effects elsewhere. Future 7-ray observations (particularly at low energies, by Fermi) could 
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help decide between these two models (PWZ and shock dominated) or explore their relative 
importance. 
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A. Appendix 

We present here the basic formulae for 77 opacity computations. Further details and 
plots (as well as the full treatment of inverse Compton interactions) can be found in the 
Appendix of Sierpowska-Bartosik & Torres (2008). The source of the thermal radiation field 
is the massive star of early type (O, Be, WR). The spectrum is described by Planck's law, 
which differential energy spectrum (the number of photons of given energy e per unit energy 
de, per unit solid angle Q, per unit volume V) is given by: 

dn(e, Si) = 4tt e 2 

dedQdV {he) 3 e^' - V 1 ' 

where e is thermal photon energy, h is the Planck constant, and k is Boltzmann constant. 

The optical depth to 7-photon absorption in the radiation field of the massive star up 
to infinity can then be calculated from the integral: 

POO 

T'yyi^Eryj X{ , OC ) / ( Ery , X{ , OC , X ry ) dX ■y , (A.2) 

Jo 

where is a photon interaction rate to e ± production in an anisotropic radiation field and 
x 7 is its propagation length. When the propagation occurs toward the massive star surface, 
the integration is performed up to the stellar surface. The photon interaction rate, A~^, is 
related to a photon of energy E 1 injected at a distance Xi from the massive star, at angle a, 
and is given by the formula: 



A 77 1 (Ey,x i ,a,Xy) = J {l + jj)dji j d(f) J ^j^^ o- 77 (/3)rfe, 



(A3) 



where x 1 is the distance to the interacting photon from the injection place along its propagat- 
ing path. The variable \i in the first integration is the cosine of the photon-photon scattering 
angle fi = cos 6. The angle <fi is the azimuthal angle between the photon (7-ray) propaga- 
tion direction and the direction to the massive star, while <p s gives its limit value for the 
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direction tangent to the massive star surface. The cross section to production is denoted 
as cr 77 (/3), where the parameter (3 in the center of mass system is (3 = \/ lo 2 — m 2 /uj, with 
u) 2 = |i? 7 e(l+cos6 l ) being the photon energy squared (in this notation it is assumed also that 
c = 1). From this we get (3 2 = 1 — 2m 2 /£ , 7 e(l + //). The kinematic condition for the angle 9 
which defines the threshold for the e ± creation is given by expression: // > nn m = ^p^- — 1. To 
simplify the equations we rewrite the internal integration in Eq. (1A3j) making use of h{fi), 
such that, 

'WiassM'k (A4) 

With the replacement (3 2 = 1 — a/e, where a = 2m 2 /E 1 (l + fi), and defining the constant 
S = 87r/(/ic) 3 , the spectrum of thermal photons is now given by the formula: 

d < e ^ = n(3) = o Q2 1 (M) 

decMdV [P> D (l-/52)2 e W(i-/3 2 )^)_i- 

Substituting in the internal integral and introducing the integration variable to ft, via 

de = 2a 0/(1 — (3 2 ) 2 d(3, yields to the integral: 

f 1 (3 a? 1 
h(a) = 2S ^ ( 1 _ (32)4e(a/{1 _ l3 2 )kTe) _ 1 ^(f3)df3. (A6) 

The lower limit of integration is from the energy condition for the process 7 + 7 —>■ e + e~ , 
i.e., it follows from the threshold condition E 7 e(l + fi) = 2m 2 , where u — m. The upper 
integration limit comes from the relativistic limit u ^> m, where we get (3 ~ 1. To proceed 
forward, we introduce a dumb variable, b, by a = kT s b, to get: 

= 2S j\kT s ) 3 b 3 ^((3) (1 _^ 2)4 e[b/{1 _] 2)) _ i d(3. (A7) 

The cross section for e ± pair production is: 

a 77 (/3) = ~r 2 7r(l - (3 2 )[(3 - (3') lni±| - 2/3(2 - (3% (A8) 

where r is the classical electron radius, and <jt — f^o * s the Thomson cross section. When 
putting this expression into the integral I\(b) (Eq. IA7I) we get finally, 



(3 _ ft i n i±| _ 2/3(2 - /3 2 ) 



(1 - /?2)3 e (6/(l-/3 2 )) _ 1 

with d = 167r(£;T s //ic) 3 (T T . 



' ' d/3, (A9) 



We parametrize the internal integral and write it as a function of b, C(b) = Ii{j3)/C\. 
Then, the integral we are after can be written as I\(b) = C\C(b) and 



In the second internal integral, we have to fix the limits (having in mind that the angle 
depends on the angle 9, then also on the variable //), 

r<t>s 



The angle 4> s determines the maximal azimuthal angle of 7-ray photon propagation with 
respect to the direction of the thermal photon, so that it gives the directions tangent to the 
star surface. This condition for <p s can be determined from an spherical triangle and the 
limits of integration with respect to the parameter /1 are the range of angles for soft photons 
coming from the star (see Sierpowska-Bartosik & Torres 2008 for details). Note that depen- 
dence on the distance to the massive star (given by the place of injection x 7 ) is already in 
the calculation of the solid angle (parameters /i and 0). When radius of the star is much 
smaller then the separation of the system d (the place of photon injection) the integral over 
/i and (ft gives the point source approximation (with the dependence (R s /d) 2 ). 
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Fig. 11. — Comparison of the optical depths in the PWZ of two different binaries: LS I +61 
303 and LS 5039. Idem as in Fig. [6] but only the model 1 was chosen for LS I +61 303 with 
the orbital solution by Grundstrom et al. (2007) and the phase of periastron arbitrary set to 
= 0. In case of LS 5039 the parameters are given in Sierpowska-Bartosik & Torres (2007). 



